subroutine srccsd(ie, fij, fab, fia, & 
          iijkl, iijka, iabij, &
          iiabj,  &
          t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, wmbej, icore, & 
          maxcore, escf, eccsd, rivv, rivo, rioo) 
use mod_iop
use mod_size
use mod_orbit
use mod_ioff
use mod_ri
implicit none
include 'omp_lib.h'
!
integer prec2, isd2, igpu2, totgpu2
common/precsoc/prec2, isd2, igpu2, totgpu2
!
integer, intent(in) :: maxcore 
real*8, intent(inout) :: escf, eccsd
real*4, intent(inout) :: ie(*), fij(*), fab(*), fia(*) 
real*4, intent(inout) :: t1(*), t2(*), r1(*), r2(*) 
real*4, intent(inout) :: iijka(*), iabij(*), iijkl(*)
real*4, intent(inout) :: fae(*), fmi(*), fme(*)
real*4, intent(inout) :: wmnij(*), wmbej(*)
real*4, intent(inout) :: iiabj(*),  icore(*)
real*4, intent(inout) :: rivv(*), rivo(*), rioo(*)

real*4, allocatable,device :: ig(:)
!
integer nsize, i, j, a, b, icrsizg
integer maxvar, maxdim, iter, gneed(5), giop(5)
integer i0, i1, i2, i3, i4, i5, i6, i7, i8, i9
integer i10, i11, i12, i13g
integer i0g, i1g, i2g, i3g, i4g, i5g, i6g, i7g, i8g, i9g, i10g, i11g, i12g
real*8 t1diag
real*8 dtmp, errdiis, eccsd2, ecc_i(200), Eccsdio(1000)
real*8 eso, tt(20), rt(20), tt1, tt2, iter_time(200)
real*4, allocatable :: idoso(:), itmp(:)
!
real*8, external :: ddot, dnrm2, nrm2so
integer, external :: idsymsoc, irpdso
!
!ldebug = .false.
!ldebug = .true.
if(totgpu==1) then 
   allocate(ig(igpusize))
endif
!
call myicopyso(iabij, r2, nioi3)
call myicopyso(fia, r1, nifvo)
!
maxvar = nifvo + nioi3 
iter = 0 
if(ldebug) then 
  call printdnrm2('iabij', iabij, 3)
endif
allocate(idoso(maxvar*10*isd), itmp(maxvar*isd))
!
!--iteration begins
!
do
iter = iter+1
if(ldebug) then 
write(6,*)
endif
write(6,*)'=============in iteration:', iter
!
call cpu_time(tt(1)); call walltimeso(rt(1))
call denomt1(t1, r1, ie, 1)
call denomt2(t2, r2, ie, 1, 1, 1)
!
maxdim = 5
!
if(ldebug) then 
  call printdnrm2('t2', t2, 3)
  write(6,*)'begin dosodiis'
endif
call cpu_time(tt1)
call dosodiis3(iter, t1, maxvar, maxdim, icore, &
              maxcore+maxvar, errdiis, idoso, itmp)
call cpu_time(tt2)
if(ldebug) then 
write(6,*)'time for dosodiis2',tt2-tt1
endif
!
call effectt2(t2, t1,  icore, 1.d0, no, nv)
!
call calcene(t1, t2, fia, iabij, eccsd)
call calcene2(t1, t2, fia, iabij, eccsd2)
write(6,*)'CCSD energy',  eccsd2
eccsdio(iter) = eccsd2
!write(6,*)'total energy', escf+eccsd, escf+eccsd2
ecc_i(iter) = eccsd2
!
call effectt2(t2, t1, icore, -1.d0, no, nv)
if(iter>maxiter) exit
!
call cpu_time(tt(2)); call walltimeso(rt(2))
write(6,'(A, F10.3, F10.3)')'====time for energy     ', tt(2)-tt(1), rt(2)-rt(1)
!===========prapere for the begin 
!iowmbej = 1
call myicopyso(fab, fae, nifvv)
call myicopyso(fij, fmi, nifoo)
call myicopyso(fia, fme, nifvo)
call myicopyso(iijkl, wmnij, nioi1)
call myicopyso(fia, r1, nifvo)
call myicopyso(iabij, r2, nioi3)
if(iowmbej==0) then 
  call myicopyso(iiabj, wmbej, nioi4)
elseif(iowmbej==1) then 
  call getallso_s(icore, 1, 'BIAJAAAS')
  call putallso_s(icore, 1, 'WMBEJAAA')

  call getallso_s(icore, 1, 'BIAJABAS')
  call putallso_s(icore, 1, 'WMBEJABA')

  call getallso_s(icore, 1, 'BIAJBBBS')
  call putallso_s(icore, 1, 'WMBEJBBB')

  call getallso_s(icore, 1, 'BIAJBABS')
  call putallso_s(icore, 1, 'WMBEJBAB')

  call getallso_s(icore, 1, 'BIAJAABS')
  call putallso_s(icore, 1, 'WMBEJAAB')

  call itranspso(icore, icore(nioi46+1), &
               vrta, popa, 1, vrtb, popb, 1, 1)
  call putallso_s(icore(nioi46+1), 1, 'WMBEJBBA')
endif
  
call cpu_time(tt(3)); call walltimeso(rt(3))
write(6,'(A, F10.3, F10.3)')'====time for prepare    ', tt(3)-tt(2), rt(3)-rt(2)
!
if(ldebug) then 
call printdnrm2('fae', fae, 11)
call printdnrm2('ijkl', iijkl, 1)
call printdnrm2('ijka', iijka, 2)
call printdnrm2('abij', iabij, 3)
call printdnrm2('rivv', rivv, 111)
call printdnrm2('rivv', rivo, 222)
call printdnrm2('rivv', rioo, 333)
call printdnrm2('iabj', iiabj, 4)
endif
giop(1:5) = 0
if(totgpu==1) then 
  call gpuneedcore(gneed, giop)
endif

!===============================================================
!==========calc for the _n5_ terms No.1  =======================
!===============================================================
!------term needed in r1 an r2 calculations 
call cpu_time(tt(4)); call walltimeso(rt(4))
if(ldebug) then 
  write(6,*)'giop', giop(1)
endif
!giop(1) = 2
if(totgpu==1.and.giop(1)==1.and.iowmbej==0) then 
   igpu=1
   i0g = 1 
   i1g = i0g + nifvo*isd 
   i2g = i1g + nioi3*isd
   i3g = i2g + nifvo*isd
   i4g = i3g + nioi3*isd
   i5g = i4g + nifvv*isd 
   i6g = i5g + nioi4*isd
   i7g = i6g + nioiri_vv*isd 
   i8g = i7g + nioiri_vo*isd
   i9g = i8g + nioiri_oo*isd
   i10g = i9g + max(nvvoo*isd*6, (nvvoo*4+nvvv*2)*isd)
   call myicopy_g(r1,    ig(i0g), nifvo)
   call myicopy_g(r2,    ig(i1g), nioi3)
   call myicopy_g(t1,    ig(i2g), nifvo)
   call myicopy_g(t2,    ig(i3g), nioi3)
   call myicopy_g(fae,   ig(i4g), nifvv)
   call myicopy_g(wmbej, ig(i5g), nioi4)
   call myicopy_g(rivv,  ig(i6g), nioiri_vv)
   call myicopy_g(rivo,  ig(i7g), nioiri_vo)
   call myicopy_g(rioo,  ig(i8g), nioiri_oo)
   call calciabci_n5_ri(ig(i0g), ig(i1g), ig(i2g),ig(i3g), &
                       ig(i4g), ig(i5g), ig(i9g), ig(i6g), ig(i7g), ig(i8g))
   call myicopy_c(ig(i0g) , r1,    nifvo)
   call myicopy_c(ig(i1g) , r2,    nioi3)
   call myicopy_c(ig(i4g) , fae,   nifvv)
   call myicopy_c(ig(i5g) , wmbej, nioi4)
!   call calciabci_n5_ri(r1, r2, t1,t2, fae, wmbej, icore, rivv, rivo, rioo)
   igpu = 0
elseif(totgpu==1.and.giop(1)==2.and.iowmbej==0) then 
   if(ldebug) write(6,*)'calciabci_n5_ri2'
   call calciabci_n5_ri2(r1, r2, t1,t2, fae, wmbej, icore, rivv, rivo, rioo, ig)
elseif((totgpu==0.or.giop(1)==0).and.iowmbej==0) then 
   if(ldebug) write(6,*)'calciabci_n5_ri'
   call calciabci_n5_ri(r1, r2, t1,t2, fae, wmbej, icore, rivv, rivo, rioo)
elseif((totgpu==0.or.giop(1)==0).and.iowmbej==1) then 
   if(ldebug) write(6,*)'calciabci_n5_ri_io'
   call calciabci_n5_ri_io(r1, r2, t1,t2, fae, icore, rivv, rivo, rioo)
elseif(totgpu==1.and.giop(1)==2.and.iowmbej==1) then 
   if(ldebug) write(6,*)'calciabci_n5_ri2_io'
   call calciabci_n5_ri2_io(r1, r2, t1,t2, fae, icore, rivv, rivo, rioo, ig)
else 
    write(6,*)'error in calc of iabci_n5'
    stop
endif

call cpu_time(tt(5)); call walltimeso(rt(5))
write(6,'(A, F10.3, F10.3)')'====time for n5 large   ', tt(5)-tt(4), rt(5)-rt(4)
if(ldebug) then 
  call printdnrm2(' fae   ', fae, 11)
  call printdnrm2(' wmbej ', wmbej, 4)
  call printdnrm2(' r1    ', r1, 22)
  call printdnrm2(' r2    ', r2, 3)
endif
call cpu_time(tt(5)); call walltimeso(rt(5))
!===============================================================
!==========calc for the abci terms No.2  =======================
!===============================================================
call cpu_time(tt(6)); call walltimeso(rt(6))
!giop(2) = 3
if(totgpu==1.and.giop(2)==1) then 
 igpu = 1
 i0g = 1 
 i1g = i0g + nifvo*isd 
 i2g = i1g + nioi3*isd
 i3g = i2g + nifvo*isd
 i4g = i3g + nioi3*isd
 i5g = i4g + nioiri_vv*isd   
 i6g = i5g + nioiri_vo*isd   
 i7g = i6g + nioiri_oo*isd 
 i8g = i7g + max(nvvoo*isd*2, nvvv*isd*2)
 call myicopy_g(r1,    ig(i0g), nifvo)
 call myicopy_g(r2,    ig(i1g), nioi3)
 call myicopy_g(t1,    ig(i2g), nifvo)
 call myicopy_g(t2,    ig(i3g), nioi3)
 call myicopy_g(rivv,  ig(i4g), nioiri_vv)
 call myicopy_g(rivo,  ig(i5g), nioiri_vo)
 call myicopy_g(rioo,  ig(i6g), nioiri_oo)
 call calciabci_ri(ig(i0g), ig(i1g), ig(i2g), ig(i3g), &
                   ig(i7g), ig(i4g), ig(i5g), ig(i6g))
 call myicopy_c(ig(i0g), r1,    nifvo)
 call myicopy_c(ig(i1g), r2,    nioi3)
 igpu = 0
elseif(totgpu==1.and.(giop(2)==2.or.giop(2)==3)) then 
 write(6,*)'giop for iabci', giop(2)
 call calciabci_ri2(r1, r2, t1, t2, icore, rivv, rivo, rioo, ig, giop)
elseif(totgpu==0.or.giop(2)==0) then 
 call calciabci_ri(r1, r2, t1, t2, icore, rivv, rivo, rioo)
else 
    write(6,*)'error in calc of iabci'
    stop
endif
call cpu_time(tt(7)); call walltimeso(rt(7))
write(6,'(A, F10.3, F10.3)')'====time for iabci      ', tt(7)-tt(6), rt(7)-rt(6)
if(ldebug) then 
  call printdnrm2(' r1    ', r1, 22)
  call printdnrm2(' r2    ', r2, 3)
endif
!testestest
!
!===============================================================
!==========calc for the abcd terms No.3 ========================
!===============================================================
if(totgpu==1.and.giop(3)==1) then 
 igpu = 1
 i0g = 1 
 i1g = i0g + nifvo*isd 
 i2g = i1g + nioi3*isd
 i3g = i2g + nifvo*isd
 i4g = i3g + nioi3*isd
 i5g = i4g + nioiri_vv*isd   
 i6g = i5g + max(nvvoo*isd*2, nvvv*isd*2)
 call myicopy_g(r1,    ig(i0g), nifvo)
 call myicopy_g(r2,    ig(i1g), nioi3)
 call myicopy_g(t1,    ig(i2g), nifvo)
 call myicopy_g(t2,    ig(i3g), nioi3)
 call myicopy_g(rivv,  ig(i4g), nioiri_vv)
 call calciabcd_ri(ig(i0g), ig(i1g), ig(i2g), ig(i3g), ig(i5g), ig(i4g))
 call myicopy_c(ig(i0g), r1,    nifvo)
 call myicopy_c(ig(i1g), r2,    nioi3)
    igpu = 0 
elseif(totgpu==1.and.(giop(3)==2.or.giop(3)==3)) then 
    call calciabcd_ri2(r1, r2, t1, t2, icore, rivv, ig, giop)
elseif(totgpu==0.or.giop(3)==0) then 
    call calciabcd_ri(r1, r2, t1, t2, icore, rivv)
else 
    write(6,*)'error in calc of iabcd'
    stop
endif
call cpu_time(tt(8)); call walltimeso(rt(8))
write(6,'(A, F10.3, F10.3)')'====time for iabcd      ', tt(8)-tt(7), rt(8)-rt(7)
if(ldebug) then 
  call printdnrm2(' r1    ', r1, 22)
  call printdnrm2(' r2    ', r2, 3)
endif
!
!===============================================================
!==========calc for the wmbej terms No.4 ========================
!===============================================================
call cpu_time(tt(9)); call walltimeso(rt(9))
!call putallso(wmbej(ioi4(1)), 1, 'WMBEJAAA')
!call putallso(wmbej(ioi4(2)), 1, 'WMBEJABA')
!call putallso(wmbej(ioi4(3)), 1, 'WMBEJBBB')
!call putallso(wmbej(ioi4(4)), 1, 'WMBEJBAB')
!call putallso(wmbej(ioi4(5)), 1, 'WMBEJAAB')
!call putallso(wmbej(ioi4(6)), 1, 'WMBEJBBA')
giop(4) = 2
if(totgpu==1.and.giop(4)==1.and.iowmbej==0) then 
 igpu = 1
 i0g = 1 
 i1g = i0g + nifvo*isd 
 i2g = i1g + nioi3*isd
 i3g = i2g + nifvo*isd
 i4g = i3g + nioi3*isd
 i5g = i4g + nioi2*isd 
 i6g = i5g + nioi3*isd
 i7g = i6g + nioi4*isd 
 i8g = i7g + max(nvvoo*isd*5, nvvv*isd*2)
 call myicopy_g(r1,    ig(i0g), nifvo)
 call myicopy_g(r2,    ig(i1g), nioi3)
 call myicopy_g(t1,    ig(i2g), nifvo)
 call myicopy_g(t2,    ig(i3g), nioi3)
 call myicopy_g(iijka, ig(i4g), nioi2)
 call myicopy_g(iabij, ig(i5g), nioi3)
 call myicopy_g(wmbej, ig(i6g), nioi4)
 call calcwmbej(ig(i4g), ig(i5g),&
                ig(i2g), ig(i3g), ig(i0g), ig(i1g), &
                ig(i6g), ig(i7g))
 call myicopy_c(ig(i0g), r1,    nifvo)
 call myicopy_c(ig(i1g), r2,    nioi3)
 igpu = 0 
elseif(totgpu==1.and.giop(4)==2.and.iowmbej==0) then 
 call calcwmbej2(iijka, iabij, &
                t1, t2, r1, r2, &
                wmbej, icore, ig)
elseif((totgpu==0.or.giop(4)==0).and.iowmbej==0) then 
 call calcwmbej(iijka, iabij, &
                t1, t2, r1, r2, &
                wmbej, icore)
elseif(iowmbej==1) then
 call calcwmbej_io(iijka, iabij, &
                t1, t2, r1, r2, &
                icore, ig)
else 
    write(6,*)'error in calc of wmbej'
    stop
endif

call cpu_time(tt(10)); call walltimeso(rt(10))
write(6,'(A, F10.3, F10.3)')'====time for wmbej      ', tt(10)-tt(9), rt(10)-rt(9)
if(ldebug) then 
 call printdnrm2(' r1    ', r1, 22)
 call printdnrm2(' r2    ', r2, 3)
endif
!===============================================================
!==========calc for the left terms No.5 ========================
!===============================================================
if(totgpu==1.and.giop(5)==1.and.iowmbej==0) then 
   igpu = 1
   i0g = 1 
   i1g = i0g +   nioi2*isd 
   i2g = i1g +   nioi3*isd
   i3g = i2g +   nioi4*isd
   i4g = i3g +   nifvo*isd
   i5g = i4g +   nioi3*isd 
   i6g = i5g +   nifvo*isd
   i7g = i6g +   nioi3*isd 
   i8g = i7g +   nifvv*isd
   i9g = i8g +   nifoo*isd
   i10g = i9g +  nifvo*isd
   i11g = i10g + nioi1*isd
   i12g = i11g + nifvo*isd
   i13g = i12g + nvvoo*4*isd
   call myicopy_g(iijka,  ig(i0g ),nioi2 )
   call myicopy_g(iabij,  ig(i1g ),nioi3 )
   call myicopy_g(iiabj,  ig(i2g ),nioi4 )
   call myicopy_g(t1,     ig(i3g ),nifvo )
   call myicopy_g(t2,     ig(i4g ),nioi3 )
   call myicopy_g(r1,     ig(i5g ),nifvo )
   call myicopy_g(r2,     ig(i6g ),nioi3 )
   call myicopy_g(fae,    ig(i7g ),nifvv )
   call myicopy_g(fmi,    ig(i8g ),nifoo )
   call myicopy_g(fme,    ig(i9g ),nifvo )
   call myicopy_g(wmnij,  ig(i10g),nioi1 )
   call myicopy_g(fia,    ig(i11g),nifvo )
   call calcleft(ig(i0g), ig(i1g), ig(i2g), ig(i3g), ig(i4g), &
               ig(i5g), ig(i6g), ig(i7g), ig(i8g), & 
               ig(i9g), ig(i10g), ig(i11g), ig(i12g)) 
   call myicopy_c(ig(i5g ),r1,     nifvo )
   call myicopy_c(ig(i6g ),r2,     nioi3 )
   igpu = 0 
elseif((totgpu==0.or.giop(5)==0).and.iowmbej==0) then 
   call calcleft(iijka, iabij, &
          iiabj, t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, fia, icore)
elseif(iowmbej==1) then 
   call calcleft_io(iijka, iabij, &
          t1, t2, r1, r2, &
          fae, fmi, fme, &
          wmnij, fia, icore)
else 
    write(6,*)'error in calc of left'
    stop
endif
call cpu_time(tt(11)); call walltimeso(rt(11))
write(6,'(A, F10.3, F10.3)')'====time for left       ', tt(11)-tt(10), rt(11)-rt(10)
if(ldebug) then 
call printdnrm2(' r1    ', r1, 22)
call printdnrm2(' r2    ', r2, 3)
endif
!if(igpuabcd==0) then 
!call myicopy_g(r1_c, icore, nifvo) 
!call axpyso(nifvo, 1.d0, icore, 1, r1, 1)
!call myicopy_g(r2_c, icore, nioi3) 
!call axpyso(nioi3, 1.d0, icore, 1, r2, 1)
!endif
!
!----
if(ldebug) then 
call printdnrm2('r1', r1, 22)
call printdnrm2('r2', r2, 3)
endif
!
write(6,'(A, F10.3, F10.3)')'====time for iteration  ', tt(11)-tt(1), rt(11)-rt(1)
iter_time(iter) = rt(11)-rt(1)
!
enddo
write(6,*)'the final energy is', Eccsd2
if(ldebug) then 
  do i = 1, iter 
     write(6,'(A, I4, 2F18.12)')'ECCSD with ITER', i, Ecc_i(i), iter_time(i)
  enddo
endif
open(555,file='CCSD_ENERGY')
write(555,*) ESCF, Eccsd2
close(555)
!
      WRITE(6,99)
99    FORMAT(T3,'     Summary of iterative solution of CC equations ')
      WRITE(6,100)
100   FORMAT(6X,59('-'))
      WRITE(6,101)
101   FORMAT(T24,'Correlation',T46,'Total',/,T8,'Iteration',T26, &
             'Energy',T46,'Energy')
      WRITE(6,100)
      DO I=1,ITER
       WRITE(6,1000) I-1, Eccsdio(I), Eccsdio(I)+ESCF, ' DIIS '
1000   FORMAT(T10,I4,T19,F18.12,T39,F19.12,T59,A6)
      ENDDO
      WRITE(6,100)

      T1DIAG = 0.D0
      T1DIAG = nrm2so(irpdso(1,nva, noa, 1), T1, 1)**2
      T1DIAG = DSQRT(T1DIAG/dble(NOA))
      WRITE(6,*) '     T1 Diagnostic =', T1DIAG
!
      WRITE(6,103) ESCF
      WRITE(6,102) ECCSD2 + ESCF
103   FORMAT('        E(SCF)            = ',F20.12)
102   FORMAT('        E(CCSD)           = ',F20.12)

return
end 
